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ABSTRACT 


This  report  summarizes  the  scientific  work  accomplished  at  the  Seismic 
Data  Analysis  Certer  (SDAC)  in  the  period  October  through  December  1974. 

During  the  quarter,  eight  technical  reports  covering  six  subjects  were 
approved  for  distribution  to  the  government  approved  list.  A letter  to  the 
editor  of  the  Bulletin  of  the  Seismological  Society  of  America  concerning 
secondary  P wave  arrivals  was  also  approved  for  publication. 

Three  of  the  reports  (TR-74-16,  17,  and  18)  describe  our  programming  and 
processing  experience  with  the  ILLIAC  IV  computing  system,  while  Report  TR-74-4 
describes  a new  source  theory  for  complex  earthquakes.  Two  reports  are 
comparative  in  nature:  the  first,  TR-74-5  compares  LASA's  Detection  Processing 

System  with  that  of  NORSAR's;  the  second,  TR-74-6  is  an  empirical  comparison 
of  the  two-segment  maximum  likelihood  frequency-wavenumber  (f-k)  spectra  with 
the  f-k  spectra  produced  by  beamforming  in  the  frequency  domain.  Studies  of 
average  P and  PKP  codas  and  of  an  iterative  beam  processor  are  documented  in 
Reports  305  and  TR-73-7,  respectively. 
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Average  P and  PKP  Coda  for  Earthquakes  (Report  305) 

Our  analysis  of  418  small-event  (n^  <.  5.8)  seismograms  recorded  at  17 
world-wide  stations,  and  of  148  large-event  (n^,  Mg  (NOS),  or  from  Pasadena 
or  Berkeley  >_  7.0)  seismograms  recorded  at  8 worldwide  stations  and  TFO 
indicates  that  coda  shape  is  primarily  a function  of  the  arrival  times  and 
relative  amplitudes  of  significant  secondary  arrivals.  However,  for  times 
greater  than  10  to  20  seconds  into  the  coda,  large-event  codas  are  approximately 
0.14  units  greater  in  amplitude  at  any  given  time  relative  to  their  maxima, 
than  the  corresponding  relative  amplitude  for  small-event  codas.  This  suggests 
that  large  events  are,  in  fact,  multiple  events,  with  the  nominal  period  of 
source  activicy  for  a given  sequence  estimated  to  be  on  the  order  of  1 to  2 
minutes.  Correspondingly,  large  events  also  tend  to  be  emergent,  displaying 
a 0.2  to  0.3  m^  increase  in  amplitude  between  5 and  30  seconds  into  the  P-wave 
arrival  over  that  observed  in  the  first  5 seconds  of  the  arrival.  Because  of 
their  differences,  large-event  and  small-event  coda  observations  cannot  be 
combined.  At  least  two  sets  of  coda  observations  are  required  (and  are 
presented  in  the  report)  for  coda  prediction.  The  small-event  cod  s are  used 
to  predict  the  codas  for  the  San  Fernando,  California,  earthquake  o. 

9 February  1971,  at  43  stations.  With  few  exceptions,  the  observed  coda  lie 
within  one  standard  deviation  of  the  predicted  coda. 


An  Iterative  Approximation  to  the  Mixed- Slgno.l  Processor  (TR-/3-7) 


In  this  study  we  develop  an  iterative  beam  processor  which  in  the  l^mit 
is  identical  to  the  mixed-signal  processor  (Dean  et.  al.,  1968).  Assuming 
that  two  events  arrive  simultaneously  at  an  array  consisting  of  N sensors, 
the  array  is  first  beamed  on  one  of  the  two  epicenters  to  produce  a signal 
estimate  for  this  event  (O’th  iteration).  This  signal  estimate  is  then  time- 
shifted  and  subtracted  om  each  of  the  original  N seismograms  in  an  attempt 
to  remove  the  signal  from  the  original  seismograms.  The  new  set  of  records, 
each  containing  N stripped  seismograms,  is  then  beamed  to  produce  a signal 
estimate  for  the  second  event.  . te  signal  estimate  for  the  second  event  is 
now  time-shifted,  subtracted  from  the  original  N se * sinograms , and  the  stripped 
seismograms  are  rebeamed  on  the  first  event.  The  process  is  repeated  until 
differences  in  successive  signal  estimates  for  the  desired  event  fall  below 
a predetermined  threshold.  The  iterative-beam  processor  has  great  practical 
'and  intuitive')  appeal.  For  seven  or  more  elements,  the  iterative  process 
converges  in  a few  iterations  requiring  only  a few  shift  and  sum  oper.  •‘iont 
per  data  point,  while  the  equivalent  mixed-signal  (asymptotic  maximum- 
likelihood)  processor  requires  a convolution  for  each  data  point. 
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A Source  Theory  tor  Complex  Earthquakes  (TR-74-4) 

!n  this  report  we  review  several  source  theories  including  that  of 
haskell  (1964)  who  developed  a theory  for  the  spectral  distribution  of  the 
teleseismic  congressional  and  shear  radiation  from  long  thin  strike-slip  or 
tonsional  earthquakes  in  a homogeneous  infinite  elastic  space.  He  noted 
that  the  strike-slip  earthquakes  had  a higher  S/P  (shear- to-compressional ) 
energy  ratio  chan  vas  commonly  observed,  and  speculated  that  most  earthquakes 
have  a small  component  of  tensional  faulting.  He  suggested  that  since 
tensional  faults  have  a lower  S/P  ratio  than  do  strike-slip  earthquakes,  a 
smaH  Portlo<i  tensional  cracking  could  dramatically  alter  the  S/P  ratio. 
Haskell  also  had  some  difficulty  in  matching  the  existing  evidence  on  the 
ratio  of  shtrt-period  to  long-period  energy,  and  speculated  that  the  fault 
-surface  might  be  "rough".  Phis  possibility  he  modeled  in  two  ways:  first, 

by  allowing  the  ramp  displacement  function  to  be  modulated  by  a sine  wave, 
and  second,  by  assuming  that  for  high  frequencies  the  fault  acts  as  the  sum 
of  several  small  liults.  Haskell’s  (1969)  pape:  developed  techniques  for 
calculating  the  nearfield  displacement  from  slip  or  tensional  earthquakes. 

In  Haskell  (1966)  he  extended  the  idea  of  a modulated  ramp  displacement 
^unction  to  a stochastic  displacement  function.  The  ensemble  of  displacement 
functions  was  characterized  by  the  average  autocorrelation  of  the  ensemble. 
Ms  average  autocorrelation  led  them  to  an  average  spectrum.  Haskell  chose 
the  autocorrelation  function  to  be  similar  to  that  of  the  sine-wave  modulated 
ramp  function  in  his  1964  paper,  except  that  it  had  no  periodic  factor,  and 
satisfied  an  integral  constraint  resulting  from  the  fact  that  the  earthquake 
must  begin  and  end  in  a state  of  static  equilibrum. 

Th?  displacement  spectrum  resulting  from  this  autocorrelation  behaved  at 
nigh  frequencies  like  3 (the  spectrum  resulting  from  the  modulated  ramp 
function  was  asymptotic  as  *“4).  Aki  ( 1967)  showed  that  an  ."3  asymptote 
leads  to  predictions  for  the  shape  of  the  curve  in  disagreement  with 

observation,  and  chose  a different  ensemble  autocorrelation  function  which 
yielded  an  asymptotic  slope  of  w-2.  (Aki’s  particular  choice  of  an  auto- 
correlation does  not  satisfy  Haskell’s  criterion  that  the  earthquake  begin 
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and  end  in  a state  of  static  equilibrium.  However,  a different  autocorrelation 
function,  which  would  satisfy  the  equilibrium  criteria,  could  be  selected  which 
would  still  result  in  an  w asymptotic  slope.) 

Savage  (1966)  extended  Haskell's  1964  nonstatistical  theory  to  handle  the 
case  of  a dislocation  nucleating  at  a point  in  the  fault  plant  and  spreading 
circularly  to  the  boundaries  of  a elliptical  fault.  (Haskell  U964)  had 
assumed  that  the  fault  started  on  a line  across  the  width  of  a narrow  fault 
and  propagated  down  the  length  of  the  fault).  Savage  showed  that  while  the 
simple  ramp  displacement  function  yielded  a displacement  spectrum  asymptotic 
as  of2  in  Haskell's  initial  line  dislocation  model;  it  resulted  in  an  u> 
asymptote  when  the  more  realistic  initial  point  dislocation  was  assumed. 

This,  he  showed,  was  due  to  the  fact  that  the  radiating  surface  grew  quadra- 
tically  with  time  instead  of  linearly.  This  «.'3  asymptote  was  corroborated 
by  Molnar,  et.  al.  (1973)  for  cases  in  which  the  fault  rupture  velocity  is 

less  than  the  shear  velocity. 

Brune  (1970)  presented  a shear  wave  source  theory  which  brought  together 
in  a physically  reasonable  way  several  important  concepts.  He  made  plausible 
the  idea  that  until  influences  from  the  ends  of  the  fault  were  felt,  the  ^ 
fault  displacement  after  cracking  would  be  a ramp  function  proportional  to  — , 
where  o is  the  stress  and  u is  the  shear  modulus.  Assuming  then,  in  effect, 
that  the  entire  fault  began  to  radiate  at  once,  the  characteristic  period  for 
the  radiated  shear  energy  became  L/B,  where  L is  the  fault  length  and  6 is  the 
sr.ear  velocity.  Brune's  shear  wave  displacement  spectrum  decayed  asymptotically 
as  u‘2.  This  slope  was  not  strictly  speaking  derived,  but  simply  followed 
from  his  assumed  form  for  the  far-field  solution  which  he  chose  to  be  in 
agreement  with  the  result  for  an  instantaneous  stress  pulse  applied  to  the 
interior  of  a spherical  surface  (see  Bullen,  1963).  This,  of  course,  is  not 
necessarily  a good  model  of  a growing  displacement  nucleating  at  a point  on 
the  surface  of  a plane  or  spherical  surface.  (The  fact  that  Bullen  [1963] 
and  Randall  [1966]  found  corner  frequencies  proportional  to  L/B  for  shear 
and  L/a  for  compressional  motion  reflects  the  fact  that  their  source  was 
assumed  to  act  instantaneously,  and  not  the  fact  that  the  source  was  distri- 
buted over  a volume  instead  of  a surface,  as  suggested  by  Molnar  et  al.  [1973]), 
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In  his  1970  paper  Brune  introduced  the  concept  of  partial  stress  drop, 
t , and  showed  that  it  could  lead  to  a spectral  slope  of  -1  somewhere  between 
zero  frequency  where  the  displacement  spectrum  is  flat,  and  asymptotic  fre- 
quencies where  an  exponent  of  -2  obtains.  Savage  (1972)  and  others  have 
pointed  out  that  for  a long  thin  fault  one  must  expect  ciW/L  where  W is  the 
fault  width  due  to  binding  of  the  fault  at  the  sides. 

Contemporaneous  with  his  theoretical  paper,  Brune  and  a number  of  co- 
workers analyzed  signals  from  many  earthquakes  and  tried  to  interpret  them 
in  its  light:  see,  for  example,  Trifunac  and  Brune  (1970),  Wyss  (1970), 

Wyss  and  Brune  (1971),  Wyss  and  Hanks  (1972),  Hanks  and  Wyss  (1972),  Molnar 
and  Wyss  (1972),  and  Wyss  and  Molnar  (1972).  One  of  the  most  interesting 
tindings  of  these  analyses  is  that  the  P wave  corner  frequency  is  higher 
than  the  S wave  corner  frequency.  Hanks  and  Wyss  (1972)  attempted  to  explain 
this  fact  by  replacing  B in  Brune's  (1970)  theory  of  a,  the  compressional 
wave  velocity,  and  appealing  to  the  concept  of  characteristic  wavelength. 

They  pointed  out,  however,  that  there  was  not  theoretical  foundation  for 
these  procedures.  Indeed  Savage  (1972)  pointed  out  that  the  more  precise 
theories  of  Haskell  (1964)  and  Savage  (1966)  predict  approximately  equal 
corner  frequencies  for  P and  S waves,  due  to  the  fact  that  the  corner  fre- 
quencies are  dominated  by  the  duration  of  faulting,  so  long  as  that  time  is 
greater  than  L/B,  and  not  by  the  length  of  the  fault  per  se. 

Molnar,  Tucker,  and  Brune  (1973)  nave  reported  that  the  P wave  corner 
frequency  is  higher  than  that  of  the  S wave  corner  for  144  aftershocks  of  the 
San  Fernando  February  9,  1971,  earthquake.  They  also  show  by  an  exact  solu- 
tion, following  Savage  (1966),  that  the  P-corner  is  higher  than  the  S-corner 
if  the  strain  is  released  instantaneously  over  the  entire  surface  of  a circular 
disk.  On  the  other  hand,  if  the  dislocation  propagates  at  about  0.5  times 
the  S wave  velocity  then  they  find  that  the  two  corner  frequencies  are  about 
equal.  This  result  is  in  agreement  with  the  studies  of  Haskell  (1964)  and 
Savage  (1966,  1972). 

At  intermediate  dislocation  velocities  the  results  of  Savage  (1974)  show 
the  importance  of  the  operational  definition  of  corner  frequency.  Savage  (1974) 
shows  that  if  the  corner  frequency  is  defined  as  the  intersection  of  the 
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high  frequency  displacement  asymptote  with  the  horizontal  low-frequency 

asymptote,  then  if  v,  the  dislocation  velocity,  is  close  to  the  shear  wave 

velocity  the  P wave  corner  will  be  substantially  higher  than  the  S wave  corner. 

On  the  other  hand,  examination  of  Savage's  Figures  1 and  2 shows  that  for  the 

same  cases  the  P wave  corner  will  be  on  the  order  of  rnly  20%  higher  than  the 

S wave  corner  if  the  corners  are  taken  as  the  6dB  (1.2  amplitude)  points  on 

the  displacement  amplitude  spectrum.  The  identical  conclusions  can  be  obtained 

by  examination  of  Figure  6 in  Molnar  et  al.  (1973).  In  fact,  one  sees  there 

the  special  dangers  of  the  asymptotic  definition  of  corner  frequencies  in 

that  for  v=a,  where  a is  the  compressional  wave  velocity,  the  intersection 

definition  gives  f <f  for  0-85°,  an  observation  point  near  the  plane  of  the 
p s 

fault.  This  result  is  easily  verified  from  Savage's  (1974)  figures,  and  one 

mav  see  also  that  for  the  6dB  corners,  f >f  . These  results  suggest  that  the 
1 p s 

6dB  corner  determination  is  more  reliable  than  determination  by  the  intersec- 
tion of  asymptotes.  It  is  probably  just  as  "objective"  as  the  intersection 
of  asymptotes  method,  and  is  certainly  not  as  severely  affected  by  Incorrect 
values  of  Q,  a serious  problem  with  the  latter  method  as  pointed  out  by 
Thatcher  and  Hanks  (1973). 

We  show  in  this  paper  that  another  explanation  of  the  f >f  paradox  may 

P s 

be  found  by  further  development  of  Haskell's  (1964,  1966,  1972)  ideas  that 

a substantial  amount  of  tensional  faulting  accompanies  many  strike-slip  or 

dip-slip  earthquakes,  and  that  many  earthquakes  are  accompanied  by  smaller 

earthquakes.  There  is  substantial  evidence  for  this  point  of  view  in  the 

literature.  Wyss  and  Brune  (1967)  concluded  that  the  28  March  1964  Alaskan 

« 

earthquake  ruptured  in  a series  of  events;  initially  propagating  in  various 
azimuthal  directions  from  the  epicenter,  then  after  a period  of  about  44 
seconds  propagating  steadily  to  the  West.  They  also  concluded  that  other 
large  earthquakes  are  similar  in  character.  Trifunac  and  Brune  (1970)  found 
that  the  Imperial  Valley  earthquakes  could  be  regarded  as  a series  of  earth- 
quakes initiated  along  a strike-slip  fault  in  accordance  with  v/B=.5-.6. 

Kamb  et  al.  (1971)  show  that  the  pattern  of  near-surface  faulting  for  the 
San  Fernando  earthquake  was  very  complicated,  and  Jung'els  and  Frazier  (1973) 
conclude  that  this  complexity  must  be  taken  into  account  in  order  to  satis- 
factorily interpret  the  static  displacement  observations.  One  might  speculate 
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that  if  the  near-surface  faulting  is  complicated,  so  too  must  be  the  deeper 
faulting.  Jungels  and  Frazier  also  require  changes  in  dip  for  the  fault 
plane.  Bolt  (1972)  analyzed  the  Pacoima  Dam  strong  motion  records  for  the 
San  Fernando  earthquake  and  concluded  that  se/eral  "bursts"  of  high  frequency 
energy  were  received  from  remote  epicenters  on  the  fault  plane.  Mikumo  (1973) 
discusses  the  substantial  field  evidence  for  flexures  in  the  San  Fernando 
fault  plane  before  adopting  a plane  surface  for  his  theoretical  model. 

Murray  (1973)  has  concluded  that  the  main  Parkfield  earthquake  consisted  of 
two  separate  earthquakes  on  two  branches  of  the  San  Andreas. 

It  seems,  in  fact,  that  almost  every  carefully  investigated  earthquake 
reveals  first-order  departures  from  a plane  slip  fault.  This,  of  course,  is 
not  surprising  given  the  inhomogeneities  of  the  real  earth.  Tensional 
sub-earthquakes  must  also  be  expected.  Although  a homogeneous  substance  may 
fail  in  shear,  one  must  expect  zones  of  weakness  in  an  inhomogeneous  material 
to  "pull  apart",  thus  introducing  a tensional  component  to  an  earthquake. 
Alternatively,  if  a fault  plane  is  curved  and  strike-slip  movement  occurs, 
one  expects  some  separation  on  the  curves. 

The  above  remarks  are  not  meant  to  deny  that  the  best  initial  approxima- 
tion to  almost  all  earthquakes  is  slip  on  a plane.  This  is  strongly  implied 
by  the  numerous  successful  studies  of  long-period  radiation  from  earthquakes 
based  on  the  double-couple  model.  However,  one  might  reasonably  expect  the 
high-frequency  displacement  spectrum  asymptotes  to  be  severely  affected  by 
small  tensional  sub-earthquakes. 

Another  puzzling  contradiction  between  theory  and  observation  has  been 
the  apparent  S/P  amplitude  ratio,  as  observed  on  velocity  or  acceleration 
seismograms,  of  about  3/1,  whereas  theory,  e.g.  Haskell  (1964),  predicts  a 
ratio  closer  to  10  or  20  to  1.  The  possibility  that  observations  are 
affected  by  unknown  variations  in  Q between  S and  P waves  somewhat  confuses 
the  issue.  The  most  reliable  data  is  presumably  that  from  closest  in,  and 
Trifunac  and  Brune  (1970)  show  strong-motion  seismograms  where  the  ratio  for 
the  1940  Imperial  Valley  sub earthquakes  vary  between  10/1  and  1/1  with  a mean 
around  3/1.  The  first  six  minutes  of  aftershocks  of  the  February  9,  1971 
San  Fernando  Earthquake  as  recorded  at  the  Pacoima  Dam  strong  motion  instrument 
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also  give  a mean  ratio  of  about  3/1,  as  we  show  in  the  original  report.  As 
pointed  out  by  Haskell  (1964),  this  paradox  could  be  resolved  if  much  of  the 
energy  around  the  corner  frequency  came  from  tensional  faulting  for  which  the 
theoretical  amplitude  ratio  is  about  3/1. 

Among  the  important  observations  which  a successful  source  theory  should 
address  are  those  which  have  grown  out  of  the  efforts  to  discriminate  earth- 
quakes from  underground  explosions.  The  amplitude  of  twenty-second  Rayleigh 
waves  plotted  versus  the  amplitude  of  one-second  comprei  sional  waves  is  known 
to  be  an  excellent  discriminant;  see  for  example  Evernden  et  al.  (1969,  1971) 

and  Marshall  and  Basham  (1972).  One  of  Aki's  (1967)  chief  motivations  in 

-3  -2 

changing  from  Haskell's  (1966)  u model  to  an  u model  was  that  the  latter 
is  in  much  tetter  agreement  with  existing  data.  On  thi  other  hand,  the 

"spread"  in  observed  M -m^  plots  is  much  greater  than  can  be  accounted  for 

by  sampling,  measurement,  and  propagation  errors;  and  therefore  any  single 

_2 

line  in  M -m^  space,  such  as  the  w model  cannot  adequately  explain  the  data. 
Examples  of  especially  discordant  earthquakes  that  look  like  explosions  with 
respect  to  this  discriminant  may  be  found  in  Landers  (1972)  and  Der  (1973). 
Douglas  et  al.  (1973a)  have  shown  that  a substantial  amount  of  this  discrepancy 
may  be  explained  by  assuming  that  the  earthquakes  are  dip-slip  with  a dip  of 
45°.  This  fault-plane  orientation  has  a radiation  pattern  maximum  for  the 
teleseismic  P radiation.  Douglas  et  al.  and  Gilvert  (1973)  have  also  shown 
that  even  for  point  sources  the  average  shallow  shearing  earthquakes  will 
separate  from  explosions  by  approximately  0. 5-1.0  magnitude  units.  A possible 
physical  explanation  is  that  if  one  first  matches  the  P wave  amplitudes,  then 
for  the  earthquake  there  is  ten  times  the  compressional  wave  amplitude  emitted 
as  shear  waves.  If  shear  waves  are  about  as  efficiently  converted  to  Rayleigh 
waves  at  the  free  surface  as  are  the  compressional  waves  and  if  the  earthquake 
is  "shallow"  with  respect  to  the  Rayleigh  wavelength,  then  one  expects  about 
ten  times  the  Rayleigh  wave  from  an  earthquake  as  from  an  explosion  of  equal 
n^.  The  same  argument  then  suggests  that  for  tensional  earthquakes  one  would 
have  only  about  three  times  the  Rayleigh  wave  amplitude. 

Anglin  (1971)  plotted  complexity  versus  third  moment  of  the  observed 
spectrum  for  earthquakes  and  explosions  and  found  good  separation.  In  general 
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there  is  no  physical  explanation  for  the  separation,  and  until  one  is  presented, 
the  weaknesses  and  limitations  of  the  discriminant  cannot  be  convincingly 
discussed.  In  this  report  we  give  a plausible  explanation  of  Anglin's  results 
by  assuming  that  the  more  complex  the  earthquake,  the  more  '-mall  earthquakes 
there  are  accompanying  it.  (A  complicating  point  in  this  discussion  is  the 
idea  developed  by  Douglf  i et  al.  (1973b)  that  complexity  may  be  due  to  arrivals 
along  separate  low  and  high-0  paths.) 

In  this  report  we  find  that  the  high-frequency  asymptotes  of  earthquakes 
-3  -1.5 

can  vary  between  w for  simple  earthquakes  to  w , where  the  exponent  is 

restrained  only  by  the  requirement  for  finite  radiated  energy  (von  Seggern 

«.nd  Blandford  [1972]  present  evidence  that  the  asymptote  for  explosions  is 
_2 

a)  ).  We  find  that  in  the  case  of  complicated  earthquakes  mixing  shear  and 

tensional  sub-earthquakes,  any  monotonically  decreasing  displacement  spectral 

shape  can  be  attained.  Savage  (1972)  has  made  much  the  same  point,  saying 

-3 

"the  spectrum  of  the  incoherent  radiation  could  mask  the  w asymptotic 
behavior".  This  implies  that  short-period  discriminants  may  be  very  unreliable 
in  practice.  However,  with  knowledge  of  the  underlying  mechanism  it  may  be 
possible  to  suitably  constrain  the  combinations  of  discriminants,  regionalize 
the  earthquakes,  and  make  short  period  discriminants  predictable  and  reliable. 

In  the  original  report,  we  develop  and  apply  the  theory,  make  some  com- 
parisons with  previous  theoretical  results,  and  compare  some  results  with 
observation. 

In  summary,  we  have  develaped  a theory  which  accounts  for  the  observed 
separation  between  Mg:m^  lines  for  earthquakes  and  explosions,  and  also 
explains  the  robustness  of  the  discriminant  despite  wide  ranges  in  earthquake 
(M  -m^)  values.  This  wide  range  is  explained  as  being  due  to  the  existence 
of  both  simple  and  complex  earthquakes.  One  also  sees  that  the  difference 
(Mg-m^)  can  be  a good  discriminant  even  though  the  slope  of  the  earthquake 
M^m^  line  is  not  1.0. 

The  fact  that  P/S  ratios  are  closer  to  0.3  than  to  0.03  as  classically 
predicted,  together  with  the  observations  of  higher  P than  S corner  frequency, 
have  been  shown  to  be  related  phenomena  and  explainable  in  terms  of  tensioral 
subearthquakes . 
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A tentative  explanation  has  been  given  for  the  success  of  the  complexity- 
third  moment  of  frequency  discriminant  (we  might  note  that  this  discriminant 
is  vulnerable  to  evasion  techniques). 

The  results  suggest  that  the  projected  intersection  of  the  earthquake 
and  explosion  populations  suggested  by  the  data  of  several  workers,  e.g. 
Marshall  and  Basham  (1972),  would  not  occur  if  data  at  lower  magnitudes  were 
obtained.  Aside  from  displacement  overshoot,  the  only  method  for  obtaining 
convergence  of  the  populations  would  be  for  small  events  to  preferentially 
draw  relaxation  energy  from  small  source  regions  as  in  the  theory  of 
Archambeau  (1968).  This  possibility  lies  outside  the  scope  of  the  present 
theory.  Mean  M^  regression  lines  appear  to  be  somewhat  lacking  in  physical 
significance.  The  fact  that  some  theoretical  justification  can  be  obtained 
for  short-period  discriminants  suggests  that  they  may  be  relied  upon  in  non- 
evasion situations  if  carefully  regionalized  and  if  care  is  taken  to  select 
paths  unaffected  by  P wave  mu. tipathing. 

Further  research  should  aim  at  deciding  if  tensional  subearthquakes 
really  exist  or  are  only  a theoretical  abstraction.  Field,  experimental,  and 
numerical  work  should  all  have  a place. 

Investigation  of  the  displacement  time  history  near  the  front  of  the 
crack  is  important  to  see  if  on  the  time  scale  of  interest  teleseismically 
(or  regionally  for  earthquake  damage  studies)  the  initial  velocity  is  dis- 
continuous as  Brune  (1970)  suggested  and  as  we  have  assumed.  Brune  (1973)  has 
presented  some  experimental  evidence  to  the  contrary,  but  it  is  not  completely 
clear  that  wave  energy  from  regions  of  earlier  movement  has  not  slightly 
distorted  the  measurements.  This  question  is,  of  course,  crucial  because  the 
discontinuity  in  the  displacement  time  history  governs  the  far-field  high- 
frequency  asymptote,  and  any  slope  substantially  different  from  uf3  would 
require  major  changes  in  the  theory. 

The  study  of  the  distribution  of  subearthquakes  is  also  crucial.  We 
have  assumed  a delta  function  here;  perhaps  a distribution  similar  to  that  for 
earthquakes  could  be  shown  to  be  reasonable  and  tested.  Again,  both  mapping 
of  faults  in  the  field  and  experimental  measurements  are  needed. 
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Although  we  do  not  expect  our  basic  results  to  change,  it  is  important 
to  replace  Haskell’s  original  solution  for  the  simple  earthquake  by  an  exact 
solution  for  a fully  two-dimensj.onal  fault  plane  surface  for  both  slip  an', 
tensional  earthquakes.  Here  Savage's  (1974)  lead  could  be  followed. 

As  a final  point,  Aki  (1967)  and  others  have  pointed  out  that  for  large 
events  ra^  is  not  a satisfactory  measure  of  the  total  1.0  Hz  energy  in  the 
signal.  A suitable  measure  suggested  by  Aki  is  the  maximum  amplitude  in  the 
P coda  times  the  "length"  of  the  coda  raised  to  some  power  between  0.5  and  1.0. 
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A Comparison  of  the  LASA  and  NORSAR  Short-Period  Arrays  (TR-74-S^ 


This  study  compares  the  LASA  and  NORSAR  short-period  arrays  in  terms  of 
their  detection  processing  systems  and  their  event  summaries,  for  data  recorded 
during  a period  of  40  days  from  15  February  to  25  March  1972.  An  overview  of 

the  worldwide  surveillance  performance  of  the  combined  LASA-NORSAR  systems  is 
also  given. 

There  are  two  signal  detection  algorithms  in  the  LASA  and  NORSAR  Detec- 
tion Processors  (DP).  The  first  algorithm  checks  successive  signal-to-noise 
ratio  threshold  crossings  by  computing  and  comparing  Short  Time  Averages  (STA) 
and  Long  Time  Averages  (LTA) . The  second  algorithm  checks  in  successive  tests 
the  consistency  of  the  azimuths  and  velocities  of  the  arriving  signal.  This 
study  showed  that  many  of  LASA/SAAC  LTA  measurements  in  the  first  algorithm 
may  include  part  of  the  signal,  thus  lowering  the  reported  signal-to-noise 
ratio.  The  LTA  measurements  in  the  NORSAR  DP  system  do  not  inlcude  any  part 
of  the  nr riving  signals. 

Either  LASA  or  NORSAR  confirmed  73%  of  the  events  on  the  NOAA  PDE 
(Preliminary  Detection  of  Epicent  .rs)  list  over  the  data  period,  and  37% 
were  confirmed  by  both  arrays.  The  LASA  alone  reported  56%  of  the  events  on 
PDE  list,  and  NORSAR  alone  reported  53%. 

V direct  comparison  of  LASA  and  NORSAR  Event  Summaries  shows  that  72%  of 
the  NORSAR  published  events  are  within  LASA's  surveillance  range.  Of  these 
in-the-range  events,  70%  were  confirmed  by  LASA.  Of  the  unconfirmed  events 
11%  were  due  to  system  failures,  and  7%  were  unconfirmed  by  DP.  Similarly, 

78%  of  the  LASA  published  events  were  within  NORSAR's  surveillance  range. 

Among  these  in-the-range  events,  38%  were  confirmed  by  NORSAR.  Of  the  uncon- 
firmed events  5%  were  due  to  system  failures,  and  45%  were  unconfirmed  by  DP. 

The  higher  percentage  of  NORSAR  DP  unconfirmed  events  is  due  partly  to  the 
high  background  noise  of  the  array. 

Although  we  do  not  estimate  the  detection  thresholds  of  the  arrays  in 
this  report,  we  note  that  the  average  noise  on  the  LASA  beams  is  about  a factor 
of  two  lower  than  that  of  NORSAR.  Therefore,  LASA's  detection  threshold  would 

be  ~0.3  magnitude  units  lower  than  NORSAR's,  if  average  signal  losses  were  the 
same  at  both  sites. 
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Comparison  of  Two  Segment  Maximum-Likelihood  (TSML)  Frequency  Wavenumber 
Spectra  with  the  Fast  Beamed  Frequency  Wavenumber  Spectra  (FKPLOT)  (TR-74-6) 

The  empirical  comparison  of  the  two  segment  maximum  likelihood  f-k  spectra 
with  the  fast  frequency  domain  beam  f-k  spectra  (FKCOMB,  FKPLOT)  shows  that 
the  latter  is  more  suitable  tor  the  separation  of  multiple  signals  it  the 
amplitudes  of  the  signals  differ  considerably.  This  advantage  of  the  beaming 
process  is  attributed  primarily  to  the  stripping  procedure,  which  has  not  been 
developed  for  the  maximum  likelihood  spectra.  The  maximum  likelihood  f-k 
spectra,  on  the  other  hand,  are  less  sensitive  to  the  array  response  and  easier 
to  interpret,  especially  if  the  array  used  contains  only  a few  elements,  since 
the  sidelobes  of  the  array  response  are  more  confusing  on  the  FKCOMB  output 
of  this  case.  Both  processes  require  a good  signal  to  noise  ratio  (-2)  for 
successful  application. 

Composited  recordings  of  long  period  Rayleigh  waves  recorded  at  LASA 
were  utilized  for  the  comparison  using  various  subarray  configurations  and 
signal  and  noise  levels. 
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A Study  of  the  ILLIAC  IV  Computer  for  Seismic  Data  Processing  (TR- 74-16) 

This  is  the  first  report  in  a series  of  three  which  pertain  to  our 
experience  with  the  ILLTAC  IV  computing  system.  The  purpose  of  this  study 
was  to  determine  the  suitability  of  the  ILLIAC  computer  for  processing  seismic 
data.  We  have  done  this  by  looking  at  the  computing  requirements  of  each  of 
several  algorithms;  and  then,  by  comparing  these  requirements  with  the 
characteristics  of  the  ILLIAC,  we  investigated  the  feasibility  of  programming 
each  of  the  algorithms  on  the  ILLIAC.  Finally,  the  procedure  FKCOMB  was 
actually  coded  for  the  ILLIAC  and  the  program  has  been  tested  and  run.  FKCOMB 
is  a long-period  seismic  signal  analysis  procedure,  which  is  important  in 
calculating  discriminants  between  earthquakes  and  nuclear  explosions;  it  may 
become  an  integral  part  of  data  processing  on  the  Seismic  Network.  FKCOMB 
was  chosen  for  this  experiment  because  the  large  amount  of  processor  time 
required  prohibits  its  use  in-house.  Also,  known  results  are  available  with 
which  to  compare  the  ILLIAC  version. 

The  ILLIAC  computer  consists  of  a control  unit,  64  arithmetic  units  or 

Q 

processing  elements  (PE),  128K  64 -bit  words  of  core,  and  10  bits  of  disk 
storage,  the  64  PE's  execute  instructions  in  lock  step;  i.e.,  they  all 
execute  the  same  instruction  simultaneously.  It  is  in  this  respect  that  the 
ILLIAC  departs  from  conventional  computer  architecture. 

The  control  unit  decodes  instructions  and  executes  instructions  for  pro- 
gram control.  It  has  24-bit  integer  arithmetic  hardware  to  calculate  indices 
and  addresses.  There  are  four  general  purpose  accumulators  and  a 64-word 
fast  access  memory  in  the  CU,  which  also  has  access  to  all  128K  of  core  and 
initiates  transfers  between  core  and  disk. 

The  primary  computational  resource  of  the  ILLIAC  is  the  array  of  64 
processing  elements.  Each  has  complete  arithmetic  capabilities  and  can 
perform  2 x 106  multiplications  per  second.  The  capacity  of  all  64  PE's  is 

g 

about  10  multiplications  per  second.  Each  PE  has  access  only  to  2K  of  core, 
and  has  only  limited  capability  to  communicate  with  other  PE's.  Control 
within  a PE  consists  of  the  ability  programatically  to  disable  selected  PE's. 
When  disabled,  a PE's  memory  is  protected  and  cannot  be  altered  by  the  PE, 
though  all  other  facets  of  instruction  are  performed. 
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The  ILL I AC  disk  is  the  primary  storage  device.  It  consists  of  13  head- 
per— track  disks  and  two  disk  controllers.  Together  the  disks  hold  approxi- 
mately 10^  bits.  Transfers  between  core  and  disk  are  initiated  by  the  ILLIAC 
control  unit  and  occur  in  blocks  or  pages  of  1024  words.  The  transfer  rate 
is  about  10^  bits  per  second.  Access  time  to  any  record  on  disk  is  40  milli- 
seconds or  less.  The  disk  can  be  loaded  from  the  Tenex  file  system  prior  to 
program  execution  and  unloaded  after  program  termination.  The  layout  of  data 
to  the  ILLIAC  disk  is  under  user  control  and  should  be  arranged  to  minimize 
access  times  Juring  program  execution. 

The  64  processing  elements  provide  the  ability  to  perform  vector  arith- 
metic operations  on  64  data  elen *nts  simultaneously.  Program  logic  generally 
requires  that  selected  PE's  be  disabled  to  avoid  redundant  calculations  if  all 
64  processing  units  are  not  required.  In  general,  program  execution  time  is 
decreased  If  disabling  of  PE's  is  avoided. 

ILLIAC  is  able  to  perform  approximately  100  million  operations  (i.e.,  a 
multiplication,  addition,  etc.)  per  second.  Any  procedure  which  requires 
fewer  than  100  billion  operations  would  have  a running  time  of  under  10  minutes. 
The  setup  time  ror  an  ILLIAC  job  is  large  enough  to  make  such  a run  impractical. 
Thus,  algorithms  requiring  very  few  computations  or  the  use  of  a small  data 
base  with  any  algorithm  are  unsuitable  for  ILLIAC  processing. 

We  conclude  from  this  study  that  the  ILLIAC  computer  when  programmed  to 
perform  seismic  processing  on  large  data  bases  can  be  a valuable  tool  in  the 
development  of  aeismic  event  detection  and  discrimination  procedures.  It  is 
feasible  to  implement  some  existing  algorithms  on  the  ILLIAC  which  are  not 
currently  used  to  process  large  data  bases,  or  some  algorithms  which  are 
proposed  but  not  tested  due  to  a lack  of  computing  power.  Our  experience 
with  one  algorithm  (FKC0MB)  which  is  representative  of  seismic  analysis  pro- 
grams shows  that  a major  benefit  of  the  ILLIAC  to  seismic  processing  is  its 
ability  to  operate  in  parallel  on  sixty-four  different  data  streams,  thereby 
reducing  the  time  required  to  process  large  data  bases.  Efficiently  arranging 
these  data  streams  for  the  processing  element  memories  is  an  important  con- 
sideration for  designing  any  seismic  algorithms  for  the  ILLIAC. 
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It  is  feasible  to  program  ILLIAC  to  perform  the  algorithms  rev  jewe*d  ir 
l.iis  study:  convolution-recursive  filtering,  PHILTRE,  matched  filtering, 

beamforming,  and  luaximum  likelihood  f-k  estimation.  Since  a uajor  factor  in 
programming  any  of  these  algorithms  is  the  data  arrangment  in  core  a ciore 
detailed  study  of  the  data  configurations  for  these  algorithms  uovi^d  be 
needed  to  optimize  the  use  of  the  computing  power  of  ILLIAC.  One  algorithm 
(FKCOMB)  was  studied  in  detail  and  implemented  on  ILIIAC  IV.  Data  editing 
schemes  were  devised  for  FKCOMB  which  can  be  used  with  appropriate  modifica- 
tions for  all  the  seismological  algorithms  we  reviewed. 

Two  independent  uses  for  ILLIAC  are  suggested.  First,  FKCOMB  and  other 
algorithms  now  used  selectively  could  be  run  routinely  on  larger  data  bases 
to  better  provide  the  services  they  already  give  on  conventional  machines. 
Second,  experimental  methods  impractical  to  test  via  conventional  machines 
could  be  tested  on  ILLIAC.  The  experience  of  implementing  FKCOMB  illustrates 
that  the  design  and  coding  of  new  algorithms  for  ILLIAC  is  not  significantly 
more  difficult  than  for  serial  machines.  The  only  phase  not  experimentally 
explored  by  this  effort  are  the  operational  problems  of  manipulating  the 
large  amounts  of  data  involved  in  routine  processing  of  long  and  short  period 
data  on  ILLIAC. 

To  maximize  efficiency,  the  time  consumed  executing  analysis  algorithms 
should  be  significantly  greater  than  the  time  required  for  data  editing.  A 
combination  of  algorithms  such  as  matched  filtering  ?nd  FKCOMB  require  an 
order  of  magnitude  more  processing  time  on  a given  memory  load  than  FKCOMB 
alone,  and  would  thereby  utilize  ILLIAC  more  efficiently. 

The  following  points  represent  our  findings  in  developing  FKCOMB  software 
on  the  ILLTAC. 

1.  Faster  and  more  reliable  network  file  transfer  between  I4-Tenex 

and  other  hosts  such  as  UCLA,  Ames,  TSS  and  SDAC  would  expedite  the  programming 
and  use  of  the  ILLIAC  system. 

2.  There  is  no  clear  preference  between  CFD  and  GLYPNIR  indicated  by 
our  experiences.  The  possibility  of  implementing  CFD  at  SDAC  or  upgrading 
service  at  UCLA  should  be  investigated,  and  an  experiment  made  in  the  use 
of  GLYPNIR  before  any  long  range  decision  is  made. 
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3.  Software  debugging  aids  presently  available  for  ILLIAC  programming  are 
minimal.  Additional  debugging  aids  would  lessen  the  task  of  ILLIAC  programming. 
User  implementation  of  such  aids  on  the  SDAC  host  is  feasible  though  at  the 
cost  of  considerable  effort. 

4.  We  estimate  that  the  time  and  effort  required  to  design  and  code 
an  ILLIAC  program  is  no  more  than  twice  that  required  for  a conventional 
machine.  Due  to  the  fact  that  ILLIAC  is  not  fully  operational  at  present 
and  the  necessity  to  handle  the  large  amounts  of  data  inherent  in  the  use 
of  an  ultrafast  machine,  the  time  and  effort  required  currently  to  debug  an 

ILLIAC  program  may  be  as  much  as  four  times  that  required  for  a conventional 
machine. 

5.  Routine  processing  of  24  hours  of  long-period  seismic  data  is  not 
feasible  at  present  due  to  the  restricted  availability  of  the  ILLIAC  processor 
and  the  inability  at  the  system  to  handle  the  large  amounts  of  data  efficiently. 


-17- 


Computer  Program  Description  (TR-74-17) 


In  this  second  ILLIAC  IV  report  we  describe  t.  preliminary  version  of  a 
long  period  array  processing  package  designed  around  the  FKCOMB  algorithm  for 
use  on  the  ILLIAC  IV  computer.  FKCOMB  is  a general-pvrpose  array-processing 
program  that  uses  frequency-wavenumber  analysis  to  pro&uce  a bulletin  which 
lists  signal  detections  and  various  statistics  for  each  detection.  Two  data 
editing  and  reformatting  modules  prepare  the  seismic  data  for  FKCOMB  and 
can  be  modified  for  use  with  other  seismic  algorithms.  Preliminary  refor- 
matting of  the  seismic  data  is  performed  by  DEMI.  The  data  is  edited  and 
fast  fourier  transformed  by  DEM2. 

The  input  parameters  required  for  operating  these  programs  and  their 
subroutines  are  described  in  this  document. 
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Computer  Listings  for  ILLIAC  IV  Version  of  FKCOMB  (TR-74-18) 


This  is  the  third  report  in  a series  of  three  published  at  the  Seismic 
Data  Analysis  Center  in  1974,  which  describe  our  studies  of  and  programming 
experience  with  the  ILLIAC  IV  computer.  The  present  report  is  a computer 
listing  of  the  ILLIAC  IV  version  of  a scientific  program  called  FKCOMB.  The 
mai  >rogum,  FKCOMB,  and  two  data-editing  and  formatting  modules,  DEMI  and 
DEM 2 were  written  in  Computational  £luid  Dynamics  Code  (CFD) ; some  subroutines 
were  written  in  ASK  code. 
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Letter  to  the  Editor  of  BSSA 


Pg  and  pP  Phases  from  Seven  Pahute  Mesa  Events 


Springer  (1974)  observed  what  he  thought  to  be  secondary  P-wave  arrivals 

in  teleseismic  records  for  Pahute  Mesa  explosions.  It  was  hypothesized  that 

these  arrivals  may  be  associated  with  a spall-closure  source  (Eisler  and 

Chilton,  1964).  Spectral  analysis  of  the  P-wave  arrivals  for  seven  Pahute 

Mesa  explosions  suggests  that  the  secondary  phase  due  to  spall  closure,  P , 

8 

may  be  present  in  the  records  for  two  of  the  events  studied. 

The  method  we  used  to  determine  the  existence  of  secondary  p arrivals 

(pP  or  P ) for  a given  event  involves  the  analysis  of  null  patterns  in  the 

P-wave  spectra  which  may  result  from  interference  of  the  types  pP-P  and  P -P. 

s 

Analysis  details  are  given  by  Cohen  (1970)  and  Cohen  et  al.  (1972). 

Analysis  of  the  records  for  seven  Pahute  Mesa  events  as  recorded  at 
five  LRSM  stations  (HN— ME,  KC-MO,  PG— BC,  RK-ON,  and  SV3QB)  and  two  observa- 
tories (CPO  and  WMO)  shows  that  the  average  P-wave  spectra  for  two  events, 
Knickerbocker  and  Rex,  exhibit  spectral-null  patterns  which  may  be  interpreted 
as  deriving  from  interference  of  the  type  P -P;  that  is,  the  null  frequencies 
fn  are  apparently  related  as  follows: 

fn  * [ (2n+l)/2]  Af,  n-0,1,2,3 

where  Af  « 1/t,  and  t is  the  P -P  delay  time.  The  spectra  for  the  five 

s 

remaining  events  display  nulls  which  appear  related  to  pP-p  interference;  here: 
fR  * nAf,  n*0,l,2,3, . . . , 

where  Af=l/x,  and  x is  the  pP-P  delay  time.  The  relationships  between  null 
freqeuncies  f^  and  integers  n for  the  seven  events  examined  are  given  in 
the  original  letter. 

Springer  (1974)  was  able  to  obtain  both  pP-P  and  P -P  delay  time  estimates 

s 

for  a given  event  because  he  used  records  taken  near  ground  surface  zero.  Our 

results,  however,  indicate  that  at  teleseismic  distances,  either  pP  or  P , 

8 

but  not  both,  may  appear  as  predominant  secondary  arij.vals.  That  the  spectra 

have  been  smoothed  by  averaging  over  a 0.15  Hz  window  cannot  account  for  our 

inability  >.0  resolve  nulls  which  may  be  associated  with  P -P  interference, 

8 
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for  even  at  delay  times  on  the  order  of  3 seconds,  the  attenuation  of 
spectral  undulations  produced  by  smoothing  is  only  about  3 db.  Rather,  it 
is  suggested  that  for  events  for  which  a spall  is  well  developed,  the  surface- 
reflected  phase  pP,  as  observed  at  teleseismic  distances,  will  be  weak.  The 
reason  for  this  is  that  if  a spall  develops,  energy  in  the  pP  phase  is  trapped 
above  the  spall,  and  is  subsequently  released  as  P . If  the  spall  Is  small 
nonexistent,  of  course,  only  pP  may  be  observed# 
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